undef("wrf_dyy_getvar")
function wrf_dyy_getvar( file_handle, varin[*]:string, time_in:integer, k_in[*]:integer, j_in[*]:integer, i_in[*]:integer )
local variable, time, var, u, v, u_in, v_in, pii, radians_per_degree, \
dims, nd, latitude, longitude, rank, kk, i, j, k

begin

i = i_in
j = j_in
k = k_in

;---As of NCL V6.0.0, wrf_dyy_get_site_uvmet can now handle a file or a list of files.
if(typeof(file_handle).eq."file") then
    ISFILE = True
    nc_file = file_handle
else
    if(typeof(file_handle).eq."list") then
        ISFILE = False
        nc_file = file_handle[0]
    else
        print("wrf_dyy_get_site_uvmet: error: the first argument must be a file or a list of files opened with addfile or addfiles")
        return
    end if
end if

variable = varin(0)
time = time_in(0)

if( (variable .eq. "uvmet") .or. (variable .eq. "uvmet10") ) then
    ;; Calculate winds rotated to earth coord. 

    pii = 3.14159265
    radians_per_degree = pii/180.

    if( (variable .eq. "uvmet") )  then
        getU = "U" 
        getV = "V" 
        if(.not. isfilevar(nc_file,"U")) then
            if(isfilevar(nc_file,"UU")) then
                getU = "UU"
                getV = "VV"
            end if
        end if
        if ( time .eq. -1 ) then
            if(ISFILE) then
                u_in = nc_file->$getU$(:,k,j,i)
                v_in = nc_file->$getV$(:,k,j,i)
            else
                u_in = file_handle[:]->$getU$(:,k,j,i)
                v_in = file_handle[:]->$getV$(:,k,j,i)
            end if
            u = wrf_user_unstagger(u_in,u_in@stagger)
            v = wrf_user_unstagger(v_in,v_in@stagger)
        else
            if(ISFILE) then
                u_in   = nc_file->$getU$(time_in,k,j,i)
                v_in   = nc_file->$getV$(time_in,k,j,i)
            else
                u_in   = file_handle[:]->$getU$(time_in,k,j,i)
                v_in   = file_handle[:]->$getV$(time_in,k,j,i)
            end if
            u = wrf_user_unstagger(u_in,u_in@stagger)
            v = wrf_user_unstagger(v_in,v_in@stagger)
        end if
    end if  ; if variable eq "uvmet"

    if( (variable .eq. "uvmet10") )  then
        if(isfilevar(nc_file,"U10")) then
            if ( time .eq. -1 ) then
                if(ISFILE) then
                    u_in   = nc_file->U10(:,j,i)
                    v_in   = nc_file->V10(:,j,i)
                else
                    u_in   = file_handle[:]->U10(:,j,i)
                    v_in   = file_handle[:]->V10(:,j,i)
                end if
                u = wrf_user_unstagger(u_in,u_in@stagger)
                v = wrf_user_unstagger(v_in,v_in@stagger)
            else
                if(ISFILE) then
                    u_in  = nc_file->U10(time_in,j,i)
                    v_in  = nc_file->V10(time_in,j,i)
                else
                    u_in  = file_handle[:]->U10(time_in,j,i)
                    v_in  = file_handle[:]->V10(time_in,j,i)
                end if
                u = wrf_user_unstagger(u_in,u_in@stagger)
                v = wrf_user_unstagger(v_in,v_in@stagger)
            end if
        else  ; may be a met file, so get lowest level of UU and VV
            if(isfilevar(nc_file,"UU")) then
                print("wrf_dyy_get_site_uvmet: Assume this is a met_em file - getting lowest level from UU and VV fields")
                if ( time .eq. -1 ) then
                    if(ISFILE) then
                        u_in  = nc_file->UU(:,0,j,i)
                        v_in  = nc_file->VV(:,0,j,i)
                    else
                        u_in  = file_handle[:]->UU(:,0,j,i)
                        v_in  = file_handle[:]->VV(:,0,j,i)
                    end if
                    u = wrf_user_unstagger(u_in,u_in@stagger)
                    v = wrf_user_unstagger(v_in,v_in@stagger)
                else
                    if(ISFILE) then
                        u_in  = nc_file->UU(time_in,0,j,i)
                        v_in  = nc_file->VV(time_in,0,j,i)
                    else
                        u_in  = file_handle[:]->UU(time_in,0,j,i)
                        v_in  = file_handle[:]->VV(time_in,0,j,i)
                    end if
                    u = wrf_user_unstagger(u_in,u_in@stagger)
                    v = wrf_user_unstagger(v_in,v_in@stagger)
                end if
            end if
        end if  
    end if ; if variable .eq. "uvmet10"


    map_projection = nc_file@MAP_PROJ


    if( any(map_projection.eq.(/0,3,6/)) ) then  ; no rotation needed
        dims = dimsizes(u)
        nd = dimsizes(dims)
        if ( nd .eq. 5 ) then
            uvmet = new( (/ 2, dims(0), dims(1), dims(2), dims(3), dims(4) /), typeof(u))
            uvmet(0,:,:,:,:,:) = u(:,:,:,:,:) 
            uvmet(1,:,:,:,:,:) = v(:,:,:,:,:)
        end if
        if ( nd .eq. 4 ) then
            uvmet = new( (/ 2, dims(0), dims(1), dims(2), dims(3) /), typeof(u))
            uvmet(0,:,:,:,:) = u(:,:,:,:) 
            uvmet(1,:,:,:,:) = v(:,:,:,:)
        end if
        if ( nd .eq. 3 ) then
            uvmet = new( (/ 2, dims(0), dims(1), dims(2) /), typeof(u))
            uvmet(0,:,:,:) = u(:,:,:) 
            uvmet(1,:,:,:) = v(:,:,:)
        end if
        if ( nd .eq. 2 ) then
            uvmet = new( (/ 2, dims(0), dims(1) /), typeof(u))
            uvmet(0,:,:) = u(:,:) 
            uvmet(1,:,:) = v(:,:)
        end if
        delete_VarAtts(u,(/"description","units"/))
        copy_VarAtts(u,uvmet)
        uvmet@description = " u,v met velocity"
        uvmet!0 = "u_v"
    end if


    if( any(map_projection.eq.(/1,2/)) ) then  ; rotation needed
        cen_lat  = nc_file@CEN_LAT
        if(isatt(nc_file,"STAND_LON")) then
            cen_long = nc_file@STAND_LON
        else
            cen_long = nc_file@CEN_LON
        end if
        true_lat1 = nc_file@TRUELAT1
        true_lat2 = nc_file@TRUELAT2

        getLAT = "XLAT" 
        getLON = "XLONG" 
        if(.not. isfilevar(nc_file,"XLAT")) then
            if(isfilevar(nc_file,"XLAT_M")) then
                getLAT = "XLAT_M"
                getLON = "XLONG_M"
            end if
        end if
        if ( time .eq. -1 ) then
            if(ISFILE) then
                latitude  = nc_file->$getLAT$(:,j,i)
                longitude = nc_file->$getLON$(:,j,i)
            else
                latitude  = file_handle[:]->$getLAT$(:,j,i)
                longitude = file_handle[:]->$getLON$(:,j,i)
            end if
        else
            if(ISFILE) then
                latitude  = nc_file->$getLAT$(time_in,j,i)
                longitude = nc_file->$getLON$(time_in,j,i)
            else
                latitude  = file_handle[:]->$getLAT$(time_in,j,i)
                longitude = file_handle[:]->$getLON$(time_in,j,i)
            end if
        end if

        cone = 1.
        if( map_projection .eq. 1) then    ; Lambert Conformal mapping
            if( (fabs(true_lat1 - true_lat2) .gt. 0.1) .and.  \
                (fabs(true_lat2 - 90. )      .gt. 0.1)       ) then
                cone = 10^(cos(true_lat1*radians_per_degree)) \
                -10^(cos(true_lat2*radians_per_degree))
                cone = cone/(10^(tan(45. -fabs(true_lat1/2.)*radians_per_degree)) - \
                10^(tan(45. -fabs(true_lat2/2.)*radians_per_degree))   )
            else
                cone = sin(fabs(true_lat1)*radians_per_degree)
            end if
        end if
        if(map_projection .eq. 2) then      ; polar stereographic
            cone = 1.
        end if
        if(map_projection .eq. 3) then      ; Mercator
            cone = 0.
        end if

        uvmet = wrf_uvmet( u, v, latitude, longitude, cen_long, cone )
        delete_VarAtts(u,(/"description","units"/))
        copy_VarAtts(u,uvmet)

    end if  

    if( (variable .eq. "uvmet10") )  then
        uvmet@description = " u10,v10 met velocity"
    end if
    return(uvmet)

end if


  if( variable .eq. "rh" ) then
       if(isfilevar(nc_file,"T")) then
         ;; first compute theta - function wrf_tk needs theta and pressure (Pa) on input
         ;; THEN compute rh, using tk, p (Pa), QVAPOR
         if ( time .eq. -1 ) then
           if(ISFILE) then
             T      = nc_file->T(:,k,j,i)
             P      = nc_file->P(:,k,j,i)
             PB     = nc_file->PB(:,k,j,i)
             QVAPOR = nc_file->QVAPOR(:,k,j,i)
           else
             T      = file_handle[:]->T(:,k,j,i)
             P      = file_handle[:]->P(:,k,j,i)
             PB     = file_handle[:]->PB(:,k,j,i)
             QVAPOR = file_handle[:]->QVAPOR(:,k,j,i)
           end if
         else
           if(ISFILE) then
             T      = nc_file->T(time_in,k,j,i)
             P      = nc_file->P(time_in,k,j,i)
             PB     = nc_file->PB(time_in,k,j,i)
             QVAPOR = nc_file->QVAPOR(time_in,k,j,i)
           else
             T      = file_handle[:]->T(time_in,k,j,i)
             P      = file_handle[:]->P(time_in,k,j,i)
             PB     = file_handle[:]->PB(time_in,k,j,i)
             QVAPOR = file_handle[:]->QVAPOR(time_in,k,j,i)
           end if
         end if  
         T = T + 300.
         P  = P + PB
         QVAPOR = QVAPOR > 0.000
         tk = wrf_tk( P , T )
         rh = wrf_rh( QVAPOR, P, tk )
         delete_VarAtts(T,(/"description","units"/))
         copy_VarAtts(T,rh)
       else
         ;; may be a met_em file - see if we can get RH
         if(isfilevar(nc_file,"RH")) then
           if ( time .eq. -1 ) then
             if(ISFILE) then
               rh = nc_file->RH(:,k,j,i)
             else
               rh = file_handle[:]->RH(:,k,j,i)
             end if
           else
             if(ISFILE) then
               rh = nc_file->RH(time_in,k,j,i)
             else
               rh = file_handle[:]->RH(time_in,k,j,i)
             end if
           end if
         end if
       end if  
       return(rh)
  end if



  if( variable .eq. "rh2" ) then
       if(isfilevar(nc_file,"T2")) then
         ;; Compute rh2, using T2, PSFC, Q2
         if ( time .eq. -1 ) then
           if(ISFILE) then
             T2   = nc_file->T2(:,j,i)
             PSFC = nc_file->PSFC(:,j,i)
             Q2   = nc_file->Q2(:,j,i)
           else
             T2   = file_handle[:]->T2(:,j,i)
             PSFC = file_handle[:]->PSFC(:,j,i)
             Q2   = file_handle[:]->Q2(:,j,i)
           end if
         else
           if(ISFILE) then
             T2   = nc_file->T2(time_in,j,i)
             PSFC = nc_file->PSFC(time_in,j,i)
             Q2   = nc_file->Q2(time_in,j,i)
           else
             T2   = file_handle[:]->T2(time_in,j,i)
             PSFC = file_handle[:]->PSFC(time_in,j,i)
             Q2   = file_handle[:]->Q2(time_in,j,i)
           end if
         end if  
         Q2 = Q2 > 0.000
         rh = wrf_rh( Q2, PSFC, T2 )
         delete_VarAtts(T2,(/"description","units"/))
         copy_VarAtts(T2,rh)
         rh@description = "2m Relative Humidity"           
       else
         ;; may be a met_em file - see if we can get RH
         if(isfilevar(nc_file,"RH")) then
           print("Probably a met_em file - get lowerst level from RH field")
           if ( time .eq. -1 ) then
             if(ISFILE) then
               rh2 = nc_file->RH(:,0,j,i)
             else
               rh2 = file_handle[:]->RH(:,0,j,i)
             end if
           else
             if(ISFILE) then
               rh2 = nc_file->RH(time_in,0,j,i)
             else
               rh2 = file_handle[:]->RH(time_in,0,j,i)
             end if
           end if
         end if
       end if  
       return(rh)
  end if


if( any( variable .eq. (/"geopt","geopotential","z","height"/) ) ) then
    ; Height [=full geopotentail height / 9.81]
    kSize = dimsizes(k)
    kk = new(kSize+1,"integer")
    kk(0:kSize-1) = k
    kk(kSize) = k(kSize-1) + 1
    if(isfilevar(nc_file,"PH")) then
        if ( time .eq. -1 ) then
            if(ISFILE) then
                var = nc_file->PH(:,kk,j,i)
                PHB = nc_file->PHB(:,kk,j,i)
            else
                var = file_handle[:]->PH(:,kk,j,i)
                PHB = file_handle[:]->PHB(:,kk,j,i)
            end if
        else
            if(ISFILE) then
                var = nc_file->PH(time,kk,j,i)
                PHB = nc_file->PHB(time,kk,j,i)
            else
                var = file_handle[:]->PH(time,kk,j,i)
                PHB = file_handle[:]->PHB(time,kk,j,i)
            end if
        end if

        var = var + PHB
        z = wrf_user_unstagger(var,var@stagger)
        z@description = "Geopotential"

    else
        ;; may be a met_em file - see if we can get GHT - data in met_em file is Height in M
        if(isfilevar(nc_file,"GHT")) then
            if ( time .eq. -1 ) then
                if(ISFILE) then
                    z = nc_file->GHT(:,kk,j,i)
                else
                    z = file_handle[:]->GHT(:,kk,j,i)
                end if
            else
                if(ISFILE) then
                    z = nc_file->GHT(time,kk,j,i)
                else
                    z = file_handle[:]->GHT(time,kk,j,i)
                end if
            end if
            z = z * 9.81
            z@description = "Geopotential"
            z@units       = "m2 s-2"
        end if
    end if

    if( any( variable .eq. (/"z","height"/) ) ) then
        z = z / 9.81
        z@description = "Height"
        z@units       = "m"
    end if
    return(z)
    delete(kk)
end if




if( any( variable .eq. (/"p","pres","pressure"/) ) ) then
    ; Full model pressure [=base pressure (PB) + pertubation pressure (P)]
    if(isfilevar(nc_file,"P")) then
        if ( time .eq. -1 ) then
            if(ISFILE) then
                var = nc_file->P(:,k,j,i)
                PB  = nc_file->PB(:,k,j,i)
            else
                var = file_handle[:]->P
                PB  = file_handle[:]->PB
            end if
        else
            if(ISFILE) then
                var = nc_file->P(time_in,k,j,i)
                PB  = nc_file->PB(time_in,k,j,i)
            else
                var = file_handle[:]->P(time_in,k,j,i)
                PB  = file_handle[:]->PB(time_in,k,j,i)
            end if
        end if
        var = var + PB
    else
        ;; may be a met_em file - see if we can get PRES
        if(isfilevar(nc_file,"PRES")) then
            if ( time .eq. -1 ) then
                if(ISFILE) then
                    var = nc_file->PRES
                else
                    var = file_handle[:]->PRES
                end if
            else
                if(ISFILE) then
                    var = nc_file->PRES(time_in,k,j,i)
                else
                    var = file_handle[:]->PRES(time_in,k,j,i)
                end if
            end if
        end if
    end if
    var@description = "Pressure"
    if( variable .eq. "pressure" ) then
        var = var * 0.01
        var@units = "hPa"
    end if
    return(var)
end if


if( any( variable .eq. (/"lat","XLAT","XLAT_M"/) ) ) then
    variable = "XLAT" 
    if(.not. isfilevar(nc_file,"XLAT")) then
        variable = "XLAT_M"
    end if
end if


if( any( variable .eq. (/"lon","long","XLONG","XLONG_M"/) ) ) then
    variable = "XLONG" 
    if(.not. isfilevar(nc_file,"XLONG")) then
        variable = "XLONG_M"
    end if
end if


if( any( variable .eq. (/"ter","HGT","HGT_M"/) ) ) then
    variable = "HGT" 
    if(.not. isfilevar(nc_file,"HGT")) then
        variable = "HGT_M"
    end if
end if


;  end of diagnostic variable list - we must want a variable already in the file.  

if ( time .eq. -1 ) then
    ;  check variable dimensionality and pull proper time  out of file
    ndims = dimsizes(filevardimsizes(nc_file,variable))
    if( ndims .eq. 4) then
        if(ISFILE) then
            var = nc_file->$variable$(:,k,j,i)
        else
            var = file_handle[:]->$variable$(:,k,j,i)
        end if
    end if
    if( ndims .eq. 3) then
        if(ISFILE) then
            var = nc_file->$variable$(:,j,i)
        else
            var = file_handle[:]->$variable$(:,j,i)
        end if
    end if
    if( ndims .eq. 2) then
        if(ISFILE) then
            var = nc_file->$variable$
        else
            var = file_handle[:]->$variable$
        end if
    end if
    if( ndims .eq. 1) then
        if(ISFILE) then
            var = nc_file->$variable$
        else
            var = file_handle[:]->$variable$
        end if
    end if
else
    ;  check variable dimensionality and pull proper time  out of file
    ndims = dimsizes(filevardimsizes(nc_file,variable))
    if( ndims .eq. 4) then
        if(ISFILE) then
            var = nc_file->$variable$(time_in,k,j,i)
        else
            var = file_handle[:]->$variable$(time_in,k,j,i)
        end if
    end if
    if( ndims .eq. 3) then
        if(ISFILE) then
            var = nc_file->$variable$(time_in,j,i)
        else
            var = file_handle[:]->$variable$(time_in,j,i)
        end if
    end if
    if( ndims .eq. 2) then
        if(ISFILE) then
            var = nc_file->$variable$(time_in,:)
        else
            var = file_handle[:]->$variable$(time_in,:)
        end if
    end if
    if( ndims .eq. 1) then
        if(ISFILE) then
            var = nc_file->$variable$(time_in)
        else
            var = file_handle[:]->$variable$(time_in)
        end if
    end if
end if

return(var)

end



undef("wrf_uv_intrp_vert")
function wrf_uv_intrp_vert(terrain:numeric, z:numeric, \
    uvmet[2][*][*]:numeric, uvmet10[2][*]:numeric, height_list[*]:numeric)
local hcol, ucol, vcol, u10, v10, alpha_u, alpha_v, uv_arr

begin


hcol = z - terrain
ucol = uvmet(0, :, :)
vcol = uvmet(1, :, :)
u10 = uvmet10(0, :)
v10 = uvmet10(1, :)

timeSize = dimsizes(u10)
heightSize = dimsizes(height_list)

; compute alpha using u10 and ucol(1)
if (any((ucol(:,1) * u10) .le. 0.0)) then
    alpha_u = 0.143 ; default value of alpha
else
    alpha_u = log(ucol(:,1)/u10(:))/log(hcol(:,1)/10.0)
end if

if ( any((vcol(:,1) * v10) .le. 0.0)) then
    alpha_v = 0.143 ; default value of alpha
else
    alpha_v = log(vcol(:,1)/v10(:))/log(hcol(:,1)/10.0)
end if

; compute alpha using ucol(0) and ucol(1)
;    if ( (ucol(1) * ucol(0)) .le. 0.0) then
;        alpha_u2 = 0.143 ; default value of alpha
;    else
;        alpha_u2 = log(ucol(1)/ucol(0))/log(hcol(1)/hcol(0))
;    end if

;    if ( (vcol(1) * vcol(0)) .le. 0.0) then
;        alpha_v2 = 0.143 ; default value of alpha
;    else
;        alpha_v2 = log(vcol(1)/vcol(0))/log(hcol(1)/hcol(0))
;    end if


; compute uv on height_list levels
;    u70 = u10 * (70.0/10.0)^alpha_u
;    v70 = v10 * (70.0/10.0)^alpha_v
;    u70_2 = ucol(0) * (70.0/hcol(0))^alpha_u2
;    v70_2 = vcol(0) * (70.0/hcol(0))^alpha_v2

u_list = new((/timeSize, heightSize/), "float")
v_list = new((/timeSize, heightSize/), "float")

do ilvl = 0, heightSize - 1
    ;   print(ilvl + " level is: " + height_list(ilvl))
    htmp = height_list(ilvl)
    if (htmp .eq. 10) then
        utmp = u10
        vtmp = v10
    else
        utmp = u10 * (htmp/10.0)^alpha_u
        vtmp = v10 * (htmp/10.0)^alpha_v
    end if
    u_list(:,ilvl) = utmp
    v_list(:,ilvl) = vtmp
end do

uv_arr = (/u_list, v_list/)
;print(uv_arr)
return(uv_arr)
end


undef("is_point_in_polygon")
function is_point_in_polygon(p[2], polygon[*][2], nCount:integer)
local nCross, x, i, p1, p2

begin
    nCross = 0
    x = new(1, "double")
    do i = 0,nCount-1,1
        p1 = polygon(i,:)
        p2 = polygon((i+1)%nCount, :)

        ; calculate the cross point of y=p.y and p1p2
        if (p1(1) .eq. p2(1)) then ; parallel to the X axis
            continue
        end if

        ; the crossover point is on the lower extension line of p1p2
        if (p(1) .lt. min((/p1(1),p2(1)/))) then
            continue
        end if
        ; the crossover point is on the lower extension line of p1p2
        if (p(1) .ge. max((/p1(1),p2(1)/))) then
            continue
        end if

        ; calculate the X of crossover point
        x = todouble(p(1) - p1(1))*todouble(p2(0)-p1(0)) / todouble(p2(1)-p1(1)) + p1(0)

        ; only count crossover point in one side
        if (x .gt. p(0)) then
            nCross = nCross + 1
        end if
    end do

    ; if nCross is even number, the point is outside of the polygon
    return ((nCross % 2) .eq. 1)
end

;===============================================================================
; function to strip the first and last quotes in strings
;
undef("strip_quotes")

function strip_quotes (instr[*]:string)

local instr, first, last, chars, outstr, strNum, i

begin
    first = 1
    strNum = dimsizes(instr)
    outstr = new(strNum, "string")
    do i = 0, strNum - 1
        chars = stringtochar (instr(i))
        last = dimsizes(chars) - 3
        outstr(i) = chartostring ( chars(first:last))
        delete(chars)
        delete(last)
    end do
    return (outstr)
end
;
; end of function strip_quotes
;===============================================================================


;===============================================================================
; function to convert irregulated time strs to udunit time
;
function timestr_to_udunits (Times:string)

local Times, split_arr

begin
    timeNum = dimsizes(Times)
    years  = new(timeNum, "integer")
    months = new(timeNum, "integer")
    days   = new(timeNum, "integer")
    hours  = new(timeNum, "integer")
    mins   = new(timeNum, "integer")
    secs   = new(timeNum, "integer")
    do i = 0, timeNum - 1
        time_str = str_strip(Times(i))
        if isStrSubset(time_str, " ") then
            split_arr = str_split(time_str, " ")
        else
            if isStrSubset(time_str, "_") then
                split_arr = str_split(time_str, "_")
            end if
        end if
        if ( dimsizes(split_arr) .eq. 1 ) then
            date_arr  = str_split(split_arr(0), "-")
            years(i)  = stringtointeger(date_arr(0))
            months(i) = stringtointeger(date_arr(1))
            days(i)   = stringtointeger(date_arr(2))
            hours(i)  = 0
            mins(i)   = 0
            secs(i)   = 0
        else
            date_arr  = str_split(split_arr(0), "-")
            time_arr  = str_split(split_arr(1), ":")
            years(i)  = stringtointeger(date_arr(0))
            months(i) = stringtointeger(date_arr(1))
            days(i)   = stringtointeger(date_arr(2))
            hours(i)  = stringtointeger(time_arr(0))
            mins(i)   = stringtointeger(time_arr(1))
            secs(i)   = stringtointeger(time_arr(2))
        end if
        delete(split_arr)
    end do
    units  = "hours since 1901-01-01 00:00:00"
    Time   = ut_inv_calendar(years,months,days,hours,mins,secs, units, 0)
    Time!0          = "Time"
    Time@long_name  = "Time"
    Time@description= "Time"
    Time@units      =  units
    Time&Time       =  Time        ; make coordinate variable
    return (Time)
end
;
; end of function timestr_to_udunits
;===============================================================================



;===============================================================================
; Interpolation of time series
;
undef("dyy_interp_ts")
function dyy_interp_ts(xi, yi, xo, opt:string)

local knots

begin
    if (opt .eq. "cubic") then ; use cubic spline
        knots = dimsizes(xi)
        return(csa1(xi, yi, knots, xo))
    end if
    if (opt .eq. "direction") then
        return(linint1(xi, yi, False, xo, 0))
    end if
    if (opt .eq. "linear") then
        return(linint1(xi, yi, False, xo, 0))
    end if
end
;
; end of function dyy_interp_ts
;===============================================================================
